B10-L10-int and B10-P10-int DE (neurons)

library(modplots)
library(tidyverse)
library(Seurat)
library(plotly)

# Mitochondrial genes
mt <- c("ENSGALG00000035334", #COX3
        "ENSGALG00000032456", #COII
        "ENSGALG00000032142", #MT-CO1
        "ENSGALG00000032079", #MT-CYB
        "ENSGALG00000037838", #ND6
        "ENSGALG00000029500", #ND5
        "ENSGALG00000036229", #MT-ND4
        "ENSGALG00000042478", #ND4L
        "ENSGALG00000030436", #ND3
        "ENSGALG00000041091", #MT-ATP6
        "ENSGALG00000032465", #MT-ATP8
        "ENSGALG00000043768", #MT-ND2
        "ENSGALG00000042750") #MT-ND1

# volcanoplot thresholds
p.adj <- 0.05
l2fc <- 0.5

B10-L10-int

# seurat object
my.se <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_lumb_int_seurat_250723.rds")
# cluster labels from B10int and L10int
ctrl_lumb_int_combined_labels <- readRDS("~/spinal_cord_paper/annotations/ctrl_lumb_int_combined_labels.rds")

identical(colnames(my.se), rownames(ctrl_lumb_int_combined_labels))

[1] TRUE

my.se$annot_sample  <- ctrl_lumb_int_combined_labels$annot_sample

my.se@active.assay <- "RNA"

intra cluster DE analysis

We do DE analysis per cluster, contrasting B10 and L10 samples:

markers <- list()
numbers <- list()
composition <-list()

for (i in seq(levels(Idents(my.se)))) {
  # subset for individual clusters
  mn.se <- subset(x = my.se, idents = levels(Idents(my.se))[i])
  mn.se$sample <- str_extract(mn.se$orig.ident, "ctrl|lumb")
  
  composition[[i]] <- mn.se[[]] %>% 
    select(sample, annot_sample)
  
  Idents(mn.se) <- "sample"
  
  tmp_markers <- FindMarkers(mn.se,
                             ident.1 = "ctrl",
                             only.pos = FALSE, 
                             min.pct = 0.25, 
                             logfc.threshold =  0.2,
                             latent.vars = c("CC.Difference.seurat"), 
                             test.use = "MAST", 
                             assay = "RNA")
  # cell numbers per sample
  numbers[[i]] <- data.frame(table(mn.se$sample))
  
  tmp_markers <- tmp_markers %>%
    rownames_to_column("Gene.stable.ID") %>% 
    left_join(gnames)
  
  markers[[i]] <- tmp_markers
}

names(markers) <- paste0("cl-", levels(Idents(my.se)))
names(numbers) <- paste0("cl-", levels(Idents(my.se)))
names(composition) <- paste0("cl-", levels(Idents(my.se)))

# bind lists into data frames
lumb_markers <- bind_rows(markers, .id = "cluster") %>% 
  mutate(cluster = factor(cluster, levels = paste0("cl-", levels(Idents(my.se)))))
lumb_numbers <- bind_rows(numbers, .id = "cluster") %>% 
  mutate(cluster = factor(cluster, levels = paste0("cl-", levels(Idents(my.se)))))
lumb_composition <- bind_rows(composition, .id = "cluster") %>% 
  mutate(cluster = factor(cluster, levels = paste0("cl-", levels(Idents(my.se)))))

saveRDS(lumb_markers, "~/spinal_cord_paper/data/Gg_ctrl_lumb_int_markers.rds")
saveRDS(lumb_numbers, "~/spinal_cord_paper/data/Gg_ctrl_lumb_int_numbers.rds")
saveRDS(lumb_composition, "~/spinal_cord_paper/data/Gg_ctrl_lumb_int_composition.rds")

# load the DE data
lumb_markers <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_lumb_int_markers.rds") %>% 
  mutate(clust_id = str_remove(cluster, "^cl-"))
lumb_numbers <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_lumb_int_numbers.rds") %>% 
  mutate(clust_id = str_remove(cluster, "^cl-"))
lumb_composition <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_lumb_int_composition.rds") %>% 
  mutate(clust_id = str_remove(cluster, "^cl-"))

Plot the cluster compositions and the number of marker genes.

DimPlot(my.se, label = TRUE, reduction = "tsne")

ggplot(data = lumb_numbers, aes(x = Var1, y = Freq, group = cluster, label = Freq)) +
  geom_col() +
  facet_wrap("cluster", nrow = 3) +
  geom_text(nudge_y = 50, size = 3) +
  ggtitle("Cluster composition by sample (nCells)")

ggplot(data = lumb_markers %>% filter(p_val_adj < 0.05 & avg_log2FC > 0.5), aes(x = cluster, group = cluster)) +
  geom_bar() +
  facet_wrap("cluster", nrow = 3,scales = "free_x") +
  ggtitle("Number of sig. DE genes")

NA
NA
NA

DA clusters from miloR

We select the clusters with the highest amount of neighbourhoods with DA of cells. Those are clusters 3 (neurons), 21, 22 (Glycogen body), and 30 motor neurons.

# vector of clusters
select_clusters <- c(3,21,22,30)

select_lumb_markers <- filter(lumb_markers, clust_id %in% select_clusters)
select_lumb_numbers <- filter(lumb_numbers, clust_id %in% select_clusters)
select_lumb_composition <- filter(lumb_composition, clust_id %in% select_clusters)

Barplots

Barplots showing number of cells from B10 or L10, and the contributions from the individual B10int and L10int clusters:

ggplot(data = select_lumb_numbers, aes(x = Var1, y = Freq, group = cluster, label = Freq)) +
  geom_col() +
  facet_wrap("cluster", nrow = 1) +
  geom_text(nudge_y = 10) +
  ggtitle("Cluster composition by sample (nCells)")

toplot <- select_lumb_composition %>% 
  mutate(annot_sample = str_extract(annot_sample, "\\d{1,2}_.{1}(?=trl|umb)")) %>% 
  mutate(annot_sample = factor(annot_sample)) %>% 
  group_by(sample, cluster) %>%
  count(annot_sample) %>% 
  mutate(label_ypos = cumsum(n)- 0.5*n)

BL_select_barplot <- ggplot(data=toplot, aes(x=sample, y=n, fill=fct_rev(annot_sample))) +
  geom_bar(stat="identity", color = "black") +
  geom_text(aes(y=label_ypos, label=annot_sample), 
            color="black", size=3.5) +
  facet_wrap("cluster", nrow = 1) +
    theme(legend.position="none") +
  ggtitle("BL10int cluster contributions")

BL_select_barplot

pdf("~/spinal_cord_paper/figures/Fig_4_BL10int_high_DA_clust_contribution.pdf")

BL_select_barplot
dev.off()

png 
  2 

Volcanoplots

toplot <- select_lumb_markers %>% 
  mutate(gene_type = case_when(
    avg_log2FC >= 0.5 & p_val_adj <= 0.05 ~ "ctrl",
    avg_log2FC <= -0.5 & p_val_adj <= 0.05 ~ "lumb",
    TRUE ~ "ns")
  )

cols <- c(ctrl = "black",
          lumb = "#419c73",
          ns = "grey")

shapes <- c(ctrl = 21,
            lumb = 21,
            ns = 20)

volplot <- ggplot(data = toplot,
                  aes(x = avg_log2FC,
                      y = -log10(p_val_adj),
                      label = Gene.name,
                      color = gene_type,
                      shape = gene_type
                  )) +
  geom_point() +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = cols) +
  scale_shape_manual(values = shapes) +
  facet_wrap("cluster", nrow = 1) +
  ylab("-log10(padj)") +
  theme_bw()+
  ggtitle("Sig. marker genes")

volplot

# ggplotly(volplot, width = 1000, height = 600)

Plots without HOX and mitochondrial genes:

# filter out HOX and MT genes
toplot_nomt <- toplot %>% 
  filter(!grepl("^HOX", Gene.name)) %>% 
  filter(!Gene.stable.ID %in% mt)

volplot_nomt <- ggplot(data = toplot_nomt,
                  aes(x = avg_log2FC,
                      y = -log10(p_val_adj),
                      label = Gene.name,
                      color = gene_type,
                      shape = gene_type
                  )) +
  geom_point() +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = cols) +
  scale_shape_manual(values = shapes) +
  facet_wrap("cluster", nrow = 1) +
  ylab("-log10(padj)") +
  theme_bw()+
  ggtitle("Sig. marker genes (without MT and HOX genes")

volplot_nomt

# ggplotly(volplot_nomt, width = 1000, height = 600)

volplot_nomt +
    ggrepel::geom_text_repel(aes(label = ifelse(gene_type == 'ns',
                                                NA,
                                                Gene.name)),
                             size = 3,
                             color = "black")

Warning: Removed 8233 rows containing missing values or values outside the scale range (`geom_text_repel()`).
Warning: ggrepel: 6 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 208 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 534 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 52 unlabeled data points (too many overlaps). Consider increasing max.overlaps

No MN DE genes

Cluster 30 (Consisting of 161 B10 vs 13 L10 cells) shows no MN related markers. Seemingly, those 13 cells are indeed motor neurons. This is supported by the fact of their expression of TUBB3, FOXP1, and SLC18A3.

# Motor neuron cluster
mn.se <- subset(x = my.se, idents = 30)
mn.se$sample <- str_extract(mn.se$orig.ident, "ctrl|lumb")

mn.se@active.assay <- "integrated"

mn_markers <- gnames[grepl("TUBB3|FOXP1|SLC18A3", gnames$Gene.name),]
mn_markers

VlnPlot(mn.se, group.by = "sample", mn_markers$Gene.stable.ID, cols = c("darkgrey", "#419c73"))

Low DA clusters from miloR

In contrast to the clusters above, we now select clusters 17, 18, and 19. The show a very low amount of DA neighbourhoods, are of similar size and represent neuronal, MFOL, and NPC populations.

# vector of MFOL clusters
low_clusters <- c(17, 18, 19)

low_lumb_markers <- filter(lumb_markers, clust_id %in% low_clusters)
low_lumb_numbers <- filter(lumb_numbers, clust_id %in% low_clusters)
low_lumb_composition <- filter(lumb_composition, clust_id %in% low_clusters)

Barplots

Barplots showing number of cells from B10 or L10, and the contributions from the individual B10int and L10int clusters:

ggplot(data = low_lumb_numbers, aes(x = Var1, y = Freq, group = cluster, label = Freq)) +
  geom_col() +
  facet_wrap("cluster", nrow = 1) +
  geom_text(nudge_y = 10) +
  ggtitle("Cluster composition by sample (nCells)")

toplot <- low_lumb_composition %>% 
  mutate(annot_sample = str_extract(annot_sample, "\\d{1,2}_.{1}(?=trl|umb)")) %>% 
  mutate(annot_sample = factor(annot_sample)) %>% 
  group_by(sample, cluster) %>%
  count(annot_sample) %>% 
  mutate(label_ypos = cumsum(n)- 0.5*n)

BL_low_barplot <- ggplot(data=toplot, aes(x=sample, y=n, fill=fct_rev(annot_sample))) +
  geom_bar(stat="identity", color = "black") +
  geom_text(aes(y=label_ypos, label=annot_sample), 
            color="black", size=3.5) +
  facet_wrap("cluster", nrow = 1) +
    theme(legend.position="none") +
  ggtitle("From which clusters do cells come?")

BL_low_barplot

pdf("~/spinal_cord_paper/figures/Fig_4_BL10int_low_DA_clust_contribution.pdf")

BL_low_barplot
dev.off()

png 
  2 

Volcanoplots

toplot <- low_lumb_markers %>% 
  mutate(gene_type = case_when(
    avg_log2FC >= 0.5 & p_val_adj <= 0.05 ~ "ctrl",
    avg_log2FC <= -0.5 & p_val_adj <= 0.05 ~ "lumb",
    TRUE ~ "ns")
  )

cols <- c(ctrl = "black",
          lumb = "#419c73",
          ns = "grey")

shapes <- c(ctrl = 21,
            lumb = 21,
            ns = 20)

volplot <- ggplot(data = toplot,
                  aes(x = avg_log2FC,
                      y = -log10(p_val_adj),
                      label = Gene.name,
                      color = gene_type,
                      shape = gene_type
                  )) +
  geom_point() +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = cols) +
  scale_shape_manual(values = shapes) +
  xlim(c(-3.5, 3.5)) +
  facet_wrap("cluster", nrow = 1) +
  ylab("-log10(padj)") +
  theme_bw()+
  ggtitle("Sig. marker genes")

volplot

# ggplotly(volplot, width = 800, height = 500)

Plots without HOX and mitochondrial genes:

# filter out HOX and MT genes
toplot_nomt <- toplot %>% 
  filter(!grepl("^HOX", Gene.name)) %>% 
  filter(!Gene.stable.ID %in% mt)

volplot_nomt <- ggplot(data = toplot_nomt,
                  aes(x = avg_log2FC,
                      y = -log10(p_val_adj),
                      label = Gene.name,
                      color = gene_type,
                      shape = gene_type
                  )) +
  geom_point() +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = cols) +
  scale_shape_manual(values = shapes) +
  xlim(c(-1.5, 1.5)) +
  facet_wrap("cluster", nrow = 1, scales = "free_y") +
  ylab("-log10(padj)") +
  theme_bw()+
  ggtitle("Sig. marker genes (without MT and HOX genes)")

volplot_nomt

# ggplotly(volplot_nomt, width = 800, height = 500)

p1 <- volplot +
    ggrepel::geom_text_repel(aes(label = ifelse(gene_type == 'ns',
                                                NA,
                                                Gene.name)),
                             size = 3,
                             color = "black")

p2 <- volplot_nomt +
    ggrepel::geom_text_repel(aes(label = ifelse(gene_type == 'ns',
                                                NA,
                                                Gene.name)),
                             size = 3,
                             color = "black")

p1

Warning: Removed 1289 rows containing missing values or values outside the scale range (`geom_text_repel()`).
Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 15 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 34 unlabeled data points (too many overlaps). Consider increasing max.overlaps

p2

Warning: Removed 1264 rows containing missing values or values outside the scale range (`geom_text_repel()`).

pdf("~/spinal_cord_paper/figures/Fig_4_BL10int_low_DA_volplots.pdf", width = 8, height = 4)

p1

Warning: Removed 1289 rows containing missing values or values outside the scale range (`geom_text_repel()`).
Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 16 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 35 unlabeled data points (too many overlaps). Consider increasing max.overlaps

p2

Warning: Removed 1264 rows containing missing values or values outside the scale range (`geom_text_repel()`).
Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider increasing max.overlaps

dev.off()

png 
  2 

# Date and time of Rendering
Sys.time()

sessionInfo()